%% Producing Figure 6: bank run dynamics

close all

load('no_policy')

%% Use cyclical components taken from HP-filter

% import data from 1980Q1 to 2023Q4, 176 quarters
% 2006Q3 is 107
data_y = readmatrix('ycig.xlsx','Sheet','Sheet1','Range','c52:fv52');
data_c = readmatrix('ycig.xlsx','Sheet','Sheet1','Range','c56:fv56');
data_i = readmatrix('ycig.xlsx','Sheet','Sheet1','Range','c64:fv64');
data_h = readmatrix('ycig.xlsx','Sheet','Sheet1','Range','c73:fv73');

log_data_y = log(data_y);
log_data_c = log(data_c);
log_data_i = log(data_i);
log_data_h = log(data_h);

hp_param = 1600;
[hp_ytrend,hp_y] = hpfilter(log_data_y,Smoothing=1600);
[hp_ctrend,hp_c] = hpfilter(log_data_c,Smoothing=1600);
[hp_itrend,hp_i] = hpfilter(log_data_i,Smoothing=1600);
[hp_htrend,hp_h] = hpfilter(log_data_h,Smoothing=1600);

total_t = size(hp_y,1);
start_t = 108;
end_t   = 108+run_before+run_after;

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','consumption: HP filter')
plot(-run_before:run_after,run2_avg_c./avg_c_nopol*100-100,'Linewidth',4)
hold on
plot(-run_before:run_after,c_b./avg_c_nopol*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,c_a./avg_c_nopol*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,hp_c(start_t:end_t)*100,'r-.','Linewidth',4)
hold on
plot(-run_before:run_after,zeros(1,1+run_before+run_after),'k:','Linewidth',2);
hold on
plot([0 0],[-30 30],'k:','Linewidth',2);
set(gca,'FontSize',40);
xlim([-run_before run_after])
xticks(-run_before:4:run_after)
ylim([-30 30]);
yticks(-30:10:30)
grid on
ytickformat('percentage')
xtickangle(0)
legend({'Model','10th-90th percentiles','','Data'},'Fontsize',30);%,'Orientation','horizontal')

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','investment: HP filter')
plot(-run_before:run_after,run2_avg_i./avg_i_nopol*100-100,'Linewidth',4)
hold on
plot(-run_before:run_after,i_b./avg_i_nopol*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,i_a./avg_i_nopol*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,hp_i(start_t:end_t)*100,'r-.','Linewidth',4)
hold on
plot(-run_before:run_after,zeros(1,1+run_before+run_after),'k:','Linewidth',2);
hold on
plot([0 0],[-30 30],'k:','Linewidth',2);
set(gca,'FontSize',40);
xlim([-run_before run_after])
xticks(-run_before:4:run_after)
ylim([-30 30]);
yticks(-30:10:30)
grid on
ytickformat('percentage')
xtickangle(0)

%% Spread, Aaa - 10-year Treasury yield and Baa - 10-year Treasury yield

% b12:b175 is 1983Q1-2023Q4.
% 2006Q3 is b106 (95).
%data_aaa    = readmatrix('AAA10Y.xlsx','Range','b12:b175');
%avg_aaa     = mean(data_aaa);
%plot_spread = data_aaa-avg_aaa;
%start_t = 95;
%end_t   = 95+run_before+run_after;

% b12:b93 is 2003Q3-2023Q4.
% 2006Q3 is b24 (13).
data_baa    = readmatrix('BAA10Y.xlsx','Range','b12:b93');
avg_baa     = mean(data_baa);
plot_spread = data_baa-avg_baa;
start_t = 14;
end_t   = 14+run_before+run_after;

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','spread')
plot(-run_before:run_after,10000*(run2_avg_spread-avg_spread_nopol),'Linewidth',4)
hold on
plot(-run_before:run_after,10000*(spread_b-avg_spread_nopol),'b:','Linewidth',3)
hold on
plot(-run_before:run_after,10000*(spread_a-avg_spread_nopol),'b:','Linewidth',3)
hold on
plot(-run_before:run_after,plot_spread(start_t:end_t)*100,'r-.','Linewidth',4)
hold on
plot(-run_before:run_after,zeros(1,1+run_before+run_after),'k:','Linewidth',2);
hold on
plot([0 0],[-200 400],'k:','Linewidth',2);
set(gca,'FontSize',40);
xlim([-run_before run_after])
xticks(-run_before:4:run_after)
ylim([-200 400]);
yticks(-200:100:400)
grid on
xtickangle(0)

%% TFP

% l239:l310 is 2006Q1-2023Q4.
% 2006Q3 is 3.
data_tfp = readmatrix('quarterly_tfp.xlsx','Sheet','quarterly','Range','l239:l310');
plot_tfp = data_tfp / 4; % translate annual rate into quarterly.
start_t = 4;
end_t   = 4+run_before+run_after;

%% TFP, take log and HP-filter

% l215:l310 is 2000Q1-2023Q4.
% 2006Q4 is 28.
% 2008Q4 is 36.
tfp_quar_chn = readmatrix('quarterly_tfp.xlsx','Sheet','quarterly','Range','l215:l310');
tfp_quar_chn = tfp_quar_chn / 400; % translate into quarterly change in natural log.
tfp_quar = ones(size(tfp_quar_chn,1),1);
for t=2:size(tfp_quar,1)
    tfp_quar(t) = tfp_quar(t-1)*exp(tfp_quar_chn(t));
end
tfp_quar = log(tfp_quar);
[tfp_trend,tfp_cyc] = hpfilter(tfp_quar,Smoothing=1600);

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','TFP')
plot(-run_before:run_after,exp(run2_avg_tfp)*100-100,'Linewidth',4)
hold on
plot(-run_before:run_after,exp(tfp_b)*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,exp(tfp_a)*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,tfp_cyc(28:28+run_before+run_after)*100,'r-.','Linewidth',4)
hold on
plot(-run_before:run_after,zeros(1,1+run_before+run_after),'k:','Linewidth',2);
hold on
plot([0 0],[-6 6],'k:','Linewidth',2);
set(gca,'FontSize',40);
xlim([-run_before run_after])
xticks(-run_before:4:run_after)
ylim([-6 6]);
yticks(-6:2:6)
grid on
ytickformat('percentage')
xtickangle(0)

%% Use CBO's potential output and hours worked

% import data from 2006Q3
% 2006Q3 is de, 2006Q4 is df
data_y_cbo = readmatrix('ycig.xlsx','Sheet','Sheet1','Range','df54:fv54');
data_c_cbo = readmatrix('ycig.xlsx','Sheet','Sheet1','Range','df62:fv62');
data_i_cbo = readmatrix('ycig.xlsx','Sheet','Sheet1','Range','df70:fv70');
data_h_cbo = readmatrix('ycig.xlsx','Sheet','Sheet1','Range','df74:fv74');

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','output: CBO potential')
plot(-run_before:run_after,run2_avg_y./avg_y_nopol*100-100,'Linewidth',4)
hold on
plot(-run_before:run_after,y_b./avg_y_nopol*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,y_a./avg_y_nopol*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,data_y_cbo(1:run_before+1+run_after)*100,'r-.','Linewidth',4)
hold on
plot(-run_before:run_after,zeros(1,1+run_before+run_after),'k:','Linewidth',3);
hold on
plot([0 0],[-30 30],'k-.','Linewidth',3);
set(gca,'FontSize',40);
xlim([-run_before run_after])
xticks(-run_before:4:run_after)
ylim([-30 30]);
yticks(-30:10:30)
grid on
ytickformat('percentage')
xtickangle(0)

figure('Units','Inches','OuterPosition',[7 0 8 8],'Name','hours worked: CBO potential')
plot(-run_before:run_after,run2_avg_h./avg_h_nopol*100-100,'Linewidth',4)
hold on
plot(-run_before:run_after,h_b./avg_h_nopol*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,h_a./avg_h_nopol*100-100,'b:','Linewidth',3)
hold on
plot(-run_before:run_after,data_h_cbo(1:run_before+1+run_after)*100,'r-.','Linewidth',4)
hold on
plot(-run_before:run_after,zeros(1,1+run_before+run_after),'k:','Linewidth',3);
hold on
plot([0 0],[-30 30],'k-.','Linewidth',3);
set(gca,'FontSize',40);
xlim([-run_before run_after])
xticks(-run_before:4:run_after)
ylim([-30 30]);
yticks(-30:10:30)
grid on
ytickformat('percentage')
xtickangle(0)

return
